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Abstract 

Recent experiments on liquid water show collective dipole orientation fluctuations dramatically 
slower then expected (with relaxation time > 50 ns) [D. P. Shelton, Phys. Rev. B 72, 020201(R) 
(2005)]. Molecular dynamics simulations of SPC/E water show large vortex-like structure of dipole 
field at ambient conditions surviving over 300 ps [J. Higo at al. PNAS, 98 5961 (2001)]. Both results 
disagree with previous results on water dipoles in similar conditions, for which autocorrelation times 
are a few ps. Motivated by these recent results, we study the water dipole reorientation using 
molecular dynamics simulations in bulk SPC/E water for temperatures ranging from ambient 300 
K down to the deep supercooled region of the phase diagram at 210 K. First, we calculate the 
dipole autocorrelation function and find that our simulations are well-described by a stretched 
exponential decay, from which we calculate the orientational autocorrelation time r a . Second, we 
define a second characteristic time, namely the time required for the randomization of molecular 
dipole orientation, the self- dipole randomization time r r , which is an upper limit on r a ; we find 
that r r ~ 5r a . Third, to check if there are correlated domains of dipoles in water which have large 
relaxation times compared to the individual dipoles, we calculate the randomization time Tb ox of 
the site-dipole field, the net dipole moment formed by a set of molecules belonging to a box of 
edge Lbox- ^ e find that the site-dipole randomization time Tb ox ~ 2.5r a for 

-^box ~ 

3A, i.e. it is 

shorter than the same quantity calculated for the self-dipole. Finally, we find that the orientational 
correlation length is short even at low T. 

PACS numbers: 
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I. INTRODUCTION 



Cooperative motion of water molecu 
both by experiments 
molecular dynamics (MD) simulations [22J, |23 




las been widely investigated in recent years, 
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When water is cooled, the cooperativity of water molecules increases. Recent experiments on 
water show large correlated domains of dipoles at ambient conditions which have a relaxation 



time much larger than the autocorrelation time of individual dipoles 



2l|. MD studies of 



water models also show the possibility of formation of large correlated domains of dipoles in 



bulk as well as interfacial water 



351 ] (where these correlated patterns of dipoles are pinned 



to solvated amino acids). These two studies are the principal motivation for the present 
investigation of the rotational cooperativity of water molecules. 

A challenging problem is to develop methods of describing molecular motion in water 
that are better able to interpret experimental results, such as incoherent quasielastic neu- 
tron scattering, light scattering, dielectric, and nuclear magnetic resonance experiments 
QQ. Severn, approximation proposals have been made for vaxious autocorrelation fnnc- 
tions describing both rotational and translational motion |20J, |27J . These methods usually 
assume the Kohlrausch- Williams- Watts stretched exponential for the long time relaxation 
behavior of autocorrelation functions <p(t), as predicted by mode coupling theory (MCT), 



4>(t) = A exp 



(1) 



The relaxation time r a , the exponent /3, and the non-ergodicit y fa ctor A are 
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param- 



eters that depend on temperature T and density p |22j, |23|, 

Our interest here is to study the orientational dynamics of water by simulating SPC/E 
water. First we calculate the orientational autocorrelation time as the fitting parameter r a 



appearing m 



Eq. © 
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Other definitions are possible, e.g., based on other fitting 



unctions for the orientational autocorrelation function decay, such as the biexponential 

in n 

26, 36J or the von Schweidler law [3J|. In all cases, the orientational autocorrelation times are 
the result of multi-parameter fitting procedures and roughly correspond to the characteristic 
time over which the orientational autocorrelation function decays by a factor of e ~ 2.7. 

To find an upper limit of the orientational autocorrelation time r a , we will introduce a 
new quantity, the dipole randomization time r r , as the time after which the fluctuations of 
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the dipoles resemble an uncorrelated random variable |3jj (Sec. IV A). We find r r > r a , and 
that r r and r a are linearly related (Sec. IV B), which is consistent with the MCT predictions 
that: 

(i) The autocorrelation times of all the autocorrelation functions of any fluctuation cou- 
pled to density fluctuations diverge at the same temperature Tmct with the same 
power law exponent; 

(ii) All the characteristic times of a supercooled liquid are proportional to one another. 

To characterize the increase of cooperativity and test for the presence of large correlated 
domains of dipoles, we also estimate the randomization time Tb ox for the site-dipole field 
(Sec. V), a quantity which measures the relaxation of the net dipole moment of all the 
molecules inside a box of edge L^, ox . Our calculations show that t^ ox when Lb ox ~ 3 A has 
a power law divergence at Tmct, but with Tb ox < T r . This result shows that the site-dipole 
field relaxes faster than the individual dipoles, resolving the apparent contradiction between 

n 

Ref. [35] and previous results. Calculations of r box for larger boxes show that r box does not 
depend on the box size and hence do not support the experimental observation of long-lived 



large domains of correlated dipoles 



21] 



II. THE SPC/E MODEL 



Our results are based on MD simulations of the extended simple point charge (SPC/E) 



model 38J]. The distance between the oxygen atom and each of the hydrogen atoms is 
0.1 nm, and the HOH angle is the tetrahedral angle 109.47° [39]. Each hydrogen atom has 
a charge qn = 0.432e, where e is the electron charge, and the oxygen atom has a charge 
qo = —2qH- I n addition, to model the van der Waals interaction, pairs of oxygen atoms of 
different molecules interact with a Lennard- Jones potential, 

12 / \ 6" 



('v../! = '< 



- - - 

Tin \Ti 



(2) 



* ,3 / \' 1,3 , 

where is the distance between molecules i and j, e = 0.65 kJ/mol and a = 0.3166 nm. 

We perform MD simulations for a system of N = 1728 molecules at density p = 1.0 g/cm 3 , 
210 K < T < 300 K, with periodic boundary conditions and a simulation time step of 1 fs. 
The temperature is controlled by the Berendsen method of rescaling the velocities 4(j] . The 



long-range Coulombic interactions |41| are treated with the reaction field technique with a 
cutoff of 0.79 nm. For each state point, we run two independent simulations to improve 
statistics. 



III. THE ORIENTATION AL AUTOCORRELATION FUNCTION d{t) 

To estimate the orientational autocorrelation time of water molecules in the supercooled 
regime, we average the scalar product of the normalized dipole vectors /Zj of each water 
molecule i in the system, 

Ci(t) = ^ &(*).&(()) 
1 - 

= ^5> os ^)>' ( 3 ) 

i=l 

where 9i(t) is the angle between jli{t) and /2i(0). This function corresponds to the average 
of the Legendre polynomial P\ (cos 9i (t)) evaluated for each molecule and can be directly 
measured by dielectric experiments. 

FigureUfa) plots C\{t) for 210 K < T < 300 K, and displays the two-step decay of typical 
glass-forming systems. The long-time regime at low T can be fit well by Eq. and the 
fitting parameters are shown in Table HI Both parameters in Eq. A and (3, show weak 
dependences on T. The resulting values of these parameters are consistent with previous 



simulations of a smaller system of SPC/E molecules 23]. 

The estimated autocorrelation times r a agree (Fig. |2J) with the power law behavior pre- 
dicted by the MCT, 

r a ~(T-T MCT )-^ . (4) 

We estimate Tmct = (194 ± 4) K and 7 a = 2.9 ± 0.3, in agreement with previous results for 
similar densities and temperatures |24j . 

The estimated values of r a , verify well the the von Schweidler law (Appendix) and the 
time-temperature superposition principle predicted by MCT, i.e. that the autocorrelation 
functions in the a-relaxation regime at different temperatures follow the same master curve 
if the time is rescaled by the autocorrelation time (Fig. [TJd) j^J. 
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IV. THE SELF-DIPOLE RANDOMIZATION TIME r r 



A. Definition and Methods 

Here we define the randomization time r r , a new quantity that we propose to characterize 
the orientational autocorrelation time. We consider the normalized dipole fa of molecule i 
over a time interval At = MSt, 

1 N 

»i = j?^2&(tk), (5) 

fe=0 

where Mi is a function of St and At, tk = kSt, and St is the time interval between two 
consecutive samples of fa. 

If St is greater than the autocorrelation time of fa, then two consecutive samples fa(t) 
and fa{t + St) are independent, hence (fa(tj) ■ faitk)) = if j ^ k, where (....) denotes the 
average over all the molecules iV in the system. Hence 

1 N 1 

= • = jfiC^ZMti) ■ = jf, (6) 



because {(fa) 2 ) = 1 for any tk, and 




fi rms = \/(tf) = -^= = \/— . (7) 



This is the result of a freely jointed chain of M bonds of the same length, for which the 
mean square end-to-end distance is N 2 (ji 2 ) = M J3]. Therefore, if St is larger than the 
orientational autocorrelation time for fa, the /i rms decreases as l/VAt. 

If, instead, St is shorter than the orientational autocorrelation time, consecutive elements 
in the sum in Eq. © are correlated (fa(t) ■ fa(t + St)) = z, resulting in a smaller fluctuation. 
This can be formally understood by considering the freely rotating chain model j2| , where 
consecutive bonds in the chain are free to rotate, each around the axis of the previous bond, 
at an angle 6, such that cos(6 l ) = z. With this assumption, the resulting mean square 
end-to-end distance for n bonds of unit length is 



(rl) = n— - 2 Zjr - fy (8) 
In the case of small 9, we have z = 1 — e + 0(e 2 ), with e = 8 2 /2 1 and z n ~ exp(— ne). 



Then, from Eq. (jHJ), we obtain 



I^Wl[2( ne -l-e-™)r • (9) 



In our problem, the bonds are dipole vectors sampled at time intervals St, and n = 
At/St = Af. Therefore Eq. (j^J) becomes 

/i rms ~^[2(AT e -l- e ^)] 1/2 . (10) 

The right-hand side of this equation behaves as l/v^At for TV — > oo, i.e., the random case 
behavior is recovered for large At /St. 

Therefore, if we define r r as the time at which the correlation goes to zero as l/v^At, it 
is possible to see that 

for any At if St > r r 

/i rms ~ 1/VAt I (11) 

I for At > r r if St < r r . 

If we consider the fluctuation of any observable, the relation (jllj) defines the randomization 
time r r for that observable |37( and r r is equal to the smallest St such that yU rms ~ 1/yAt 
for any At. 



B. Calculation of r r 

In Fig. 01 we show yU rms for T = 220 K calculated for different values of St. For small 
St and small At, /i rms deviates largely from the asymptotic law. However, for increasing St, 
the deviation decreases. For St = 1280 ps the asymptotic behavior, within the error of our 
calculations, is reached. 

The evaluation of r r from a plot such as in Fig. El could be problematic, since it depends 
critically on the data errors. Therefore, to define in a clear way r r , we fit the first eight 
points (At = St, 2St, 8St) using 

/i rms ~ (At)\ (12) 

where A = X(St). In this way we study how the deviation from the asymptotic regime 
decreases by increasing St. We find that the exponent A increases toward the asymptotic 
value 1/2 for increasing St, and A = 1/2 for any St > r r (Fig. 0]). We therefore define r r as 
the extrapolated values of St at which A = 1/2. We find that A approaches 1/2 as 1/St, to 
the leading order, for low temperatures (Fig. EJ). 

The resulting values of r r are presented in Fig|I)k as functions of T — T MCT , showing 
that the power law behavior Eq. (HJ) is well satisfied by r r . In this case our estimates are 
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Tmct = (191.5 ± 2.5) K and 7 r = 3.3 ± 0.2, both consistent within the errors with the 
estimates based on r a (Fig|2J). Therefore, the prediction (i) of MCT is verified. 

By plotting r r against T a , we verify the MCT prediction (ii). We find (Fig. EJ) that r r and 
r a are linearly related and that r T is approximately five times larger than r a . 

The large value of r r with respect to r a is consistent with the fact that the latter measures 
the decay of the self-dipole correlation to a finite value, while the first measures the time 
needed for the self-dipole autocorrelation to decay to zero. This result is also reminiscent 
of the recent MD analysis in bulk water for the site-dipole field, a measure of the average 



orientation of the mo 
Ref. 



351 ] . Higo et al. 



ecules passing through each spatial position, recently introduced in 
find coherent patterns for the site-dipole field, at ambient pressure 
and T = 298 K and T = 300 K, that persist for more than 100 ps, a time much larger than 
the single molecule orientational relaxation time r a of approximately 5 ps (Table I). It is, 
therefore, interesting to calculate the randomization time Tb ox and to find its relation with 
the autocorrelation time r a for T — > Tmct- 

V. THE SITE-DIPOLE FIELD 

To check if there are large correlated domains of dipoles in water which have large re- 
laxation times compared to the individual dipole correlation time, we next study site-dipole 
field introduced by Higo et. al. We define the instantaneous coarse-grained site-dipole 
field 

^ > box 

as the average of dipoles jli of all the molecules rii(t) at time t belonging to box i of edge 
£box? volume v = L^ ox and centered at fl. If rii(t) = 0, then d\ = by definitio n |43l l44j . 
We chose vectors f\ in such a way that the corresponding boxes do not overlap |45j. The 
time average d\ over an interval At, is defined analogously to Eq. (jSJ). The rms average firms; 
is defined in analogy to Eq. (0) and (JTj), but instead of summation over all molecules we 
perform a summation over all boxes. 

Since the argument presented for /i rms is also valid for c^m S > the relation (fTT|) also holds 
for d" ms and allows us to estimate the randomization time r box for d v vms . We find that r box , 
calculated for Lb ox = 3.33A , diverges at T^ ox = (194 ±2) K with a power law with exponent 
Ibox = 3.2 ± 0.2, consistent with our estimates of 7 a and Tmct, respectively (Fig. |7j). 
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If we compare Tb ox with r r (Fig. |HJ), we again find a linear relation, as in Fig. El for r a , 
consistent with the MCT statement (ii). The proportionality factor is approximately 2.5 



46| , smaller than the factor « 5 found for r r in Fig. El Therefore, we conclude that in bulk 



water the coarse-grained site-dipole randomization time Tb ox is larger than the self-dipole 
autocorrelation time r a , but smaller than r r . Thus we do not find a significant increase in 
the box dipole autocorrelation time compared to the autocorrelation time r a . These results 
do not support the results of Refs. (21I [3^ . 

To test the existence of cooperative domains in the SPC/E model, we perform coarse- 
graining of the dipole field for boxes of sizes 3.33 A < Lb ox < 10 A. If the dipoles of 
molecules in the box are independent random variables, d^ ms must be inversely proportional 
to y/< Wi > oc \/v, since the average number of molecules in the box is proportional to its 
volume. The dependence of d v rras y/v on time t must be the same for the boxes of different 
volumes v. We show in Fig. El the behavior of d^^^/v for T = 220 K and T = 300 K. 
The collapse of all the curves confirms the hypothesis of very weak autocorrelations among 
neighboring dipoles. Only for T = 220 K we observe a weak size dependence of d^ ms y/v for 
the smallest size, suggesting that at this T the correlation length is between 3.33 and 6 A, 
comparable to the dipole-dipole correlation length at ambient T 44]. Thus our simulations 
support the existence of only short range orientational autocorrelation in SPC/E water even 
at low T. 



VI. DISCUSSION 

Considerable numerical evidence shows that MCT predictions apply to orientational dy- 
namics of water, despite the fact that MCT has been developed for particles interacting 
through spherically symmetric potentials |4j|. However, recent extensions of MCT to lie 




uids of linear molecules (4$, |49j , and single solute molecules in a simple solvent liquic 
confirm the main MCT predictions about the orientational autocorrelation functions 

Our study of supercooled water confirms the validity of MCT predictions for the orienta- 
tional autocorrelation time r a , estimated through a stretched exponential of the dipole auto- 
correlation function, for the temperature range 210 K< T < 300 K at density p = 1 g/cm 3 . 
Our results agree with the time-temperature superposition principle and the power law 
Eq. ©, with T MC T = (194 ± 4) K and 7 a = 2.9 ± 0.3. 
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By evaluating the randomization time r r , defined as the time needed to randomize the 
molecular dipoles, we verify the MCT prediction that all the characteristic times of quantities 
coupled to density fluctuations of a supercooled liquid are proportional to each other and 
follow the same power law Eq. We find r r w 5 r tt , with Tmct = (191.5 ± 2.5) K and 
7 r = 3.3 ± 0.2, consistent with the estimates based on the calculation of r a . 

We also calculate the randomization time Tb ox for the box dipole field, a quantity intro- 
duced in Ref. 3^| to measure the local orientational memory of molecules passing through a 
given spatial position. Our results for L box = 3.33 Ashow that r box diverges at T box = T MCT , 
following a power law with exponent 7b ox = 7a, and that Tb ox ~ r r /2. As a consequence, the 
local memory is lost faster than the self-dipole orientational memory. 

Our results also show the existence of domains of correlated dipoles of short spatial range, 
with a correlation length comparable to the dipole-dipole correlation length at ambient T 
Q|. Whether this conclusion is specific to the SPC/E model with reaction field is an open 
question, and requires further investigation using other models of water, e.g. polarizable 
models. 
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APPENDIX A: THE VON SCHWEIDLER LAW 



The MCT predicts that the autocorrelation function departs from the plateau A as a 
power law with exponent b, known as the von Schweidler law, 

Cl (t)-A~-(t/T a f , (Al) 

where the von Schweidler exponent h does not depend on T. We verify that at lower 
temperatures Eq. (jAljl holds for roughly two decades in time (Fig. ITUj) and we find a clear 
deviation only for T > 260 K at short times, possibly due to the fact that for T > 260 K it 
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is more difficult to estimate the plateau A. The estimated value of b is 0.6 ± 0.1, consistent 
with previous results j^j and with the MCT prediction that 7 a , a, and b are related by the 
equation 

la = — + 4- (A2) 

la 2a 2b K J 

Here a is the exponent of the power law that describes the short-time approach to the plateau 
C\ — A ~ t~ a , and a is related to b by the transcendental equation 

[r(i - a)] 2 _ [r(i + b)f 

T(l-2a) r(l + 26) ' 1 ' 

where T(x) is the Euler gamma function. Our estimates of b and 7 a are consistent with both 
Eqs. (JMJ) and (ESJ) with a = 0.25 ± 0.05. 

The values of the exponents a, b and 7 a are not universal, but depend on density. However, 
the rescaling of the autocorrelation functions for different T on the same master curve, shows 
that the orientational correlation function depends on T and p only through the dependence 
on r a , as predicted by the MCT. 
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TABLE I: Parameters of the fit of C\{t) in Fig. ^ w ith Eq. (^Q). The error on each parameter is 

±10%. 

T(K) A r^(ps) /3_ 

300 0.93 4.9 x 10° 0.88 

260 0.94 1.7 x 10 1 0.85 

250 0.94 2.8 x 10 1 0.85 

240 0.94 4.9 x 10 1 0.84 

230 0.94 1.1 x 10 2 0.84 

220 0.94 2.7 x 10 2 0.83 

210 0.94 1.1 x 10 3 0.82 
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FIG. 1: (a) The orientational autocorrelation function C\ as a function of time t for T(K) = 210 
(O), 220 (□), 230 (0), 240 (A), 250 (<), 260 (y), 300(>). Symbols are simulations, lines are fits 
over the range for t > 0.03 ps to Eq. (^Q) with the fitting parameters listed in Table [IJ (b) Test of 
the time-temperature superposition principle, as predicted by MCT. The symbols and the lines for 
different T fall on a single curve if the times are rescaled by r a (T). 
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FIG. 2: (a) Power law behavior of orientational autocorrelation time r a extracted from Ci(t), as 
a function of T — Tmct> as predicted by MCT [Eq. The line is a fit to the MCT power law 
with Tmct = 194 K and exponent 7 a = 2.9. (b) To optimize the fit, we vary Tmct, calculate 
the autocorrelation coefficient R (solid line) and the x 2 (dashed line), and choose as our estimate 
of Tmct the value corresponding to the maximum or minimum of these quantities, within a 10% 
variation in our range of T. R and x 2 are rescaled to the maximum and minimum values we found 
for 188 K < Tmct < 202 K. (c) The MCT exponent j a corresponding to different choices of Tmct- 
Note that the exponent "f a decreases almost linearly with increasing choice of Tmct- Based on the 
results in (b), our estimates are Tmct = (194 ± 4) K and 7 a = 2.9 ± 0.3. 
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FIG. 3: (a) The 

Mrms of Eq. for T — 220 K, plotted versus At for a range of different time 
steps St = 128 ps (O), 576 ps (+), 832 ps (x), 1088 ps (A), and 1280 ps (□). Dashed lines show 
the predicted asymptotic behavior /i rms ~ l/y/At. The fit with Eq. (fTU|) (solid lines) is good when 
St > 576 ps, but we are unable to fit the data for 5t = 128 ps, showing that the angle between 
the dipoles in Eq. © is not independent, as assumed in the freely rotating chain model. However, 
Eq. ()10|) gives a fair description of the approach to the asymptotic regime, (b) n rms \^At vs At 
approaches a constant asymptotically when At 3> r r . In both panels (a) and (b) the errors are 
roughly the size of the symbols. 
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FIG. 4: The exponent A, denned in Eq. (|12|) and calculated using the first eight points of the curves 
in Fig. 01 versus the inverse of time-step St, for T = 300 K (O), 260 K (■), 250 K (0), 240 K 
(a), 230 K (<), 220 K (▼). Where not shown the errors are smaller than the symbol size. The 
horizontal dashed line corresponds to A = 0.5. By a quadratic fit of the data with A < 0.5, we find 
the self-dipole randomization time r r , defined as the value of 5t where A = 1/2. 
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FIG. 5: Analog of Fig. 2 for the self-dipole randomization time r r . (a) r r follows a power law 
behavior in T — Tmct [Eq. ©], as predicted by MCT. The line is a fit with Tmct = 191.5 K and 
power 7 r = 3.3. (b) As in Fig. El to optimize the estimate of Tmct we calculate the autocorrelation 
coefficient R (solid line) and the x 2 (dashed line). In (b), R and x 2 are rescaled to the maximum and 
minimum values we found for 187 K < Tmct < 200 K. (c) The fitting parameter 7 r corresponding 
to different estimates of Tmct- The exponent j r decreases linearly with increasing estimates of 
Tmct- Based on the results in (b), our estimates are Tmct = (191.5 ± 2.5) K and j r = 3.3 ± 0.2. 
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FIG. 6: Parametric plot of the times r r (T) and r (T), within the range 220 K < T < 300 K, 
with the lowest time corresponding to the highest T. The line reflects the linear one-parameter fit, 
r r = (5.1±0.2)r o . 
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FIG. 7: Analog of Fig. 2 and Fig. 6 for the site-dipole randomization time Tb ox . (a) We find a 
power law behavior in T — Tb ox , calculated for L^ ox = 3A. The line is a fit with Tmct = 194 K and 
exponent 7b ox = 3.2. (b) Optimization analysis for Tb ox : correlation coefficient R (solid line) and \ 2 
(dashed line), both rescaled to the maximum and minimum values found for 188 K< Tb ox < 202 K. 
(c) The exponent 7b ox corresponding to different choices of Tb ox , decreases linearly with increasing 
choice of T box . We estimate T box = (194 ± 2) K and 7 box = 3.2 ± 0.2. 
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FIG. 8: Analog of Fig. 7: A parametric plot of the site-dipole randomization time t^ ox (T) and the 
orientational autocorrelation time T a (T) over the range 220 K < T < 300 K, with the lowest time 
corresponding to the highest T. The line reflects the linear one-parameter fit Tb ox = (2.5 ± 0.2)r„. 
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FIG. 9: Size dependence of site-dipole autocorrelation function d^ ms as a function of t for Lbox = 
3.33 A (O), 6 A (□), 10 A (0) and for two different temperatures (a) T = 300 K and (b) 
T = 220 K. In (a) the line is a fit of data for L box = 3.33 A with < ms = a/Vi, with a = 0.08 ±0.01. 
In (b) the same fit is for the data at Lbox = 6 A and t > 10 2 , with a = 0.49 ± 0.01. For each T 
all the values of drmsV^ overlap, suggesting that the orientational autocorrelation is short-range. 
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FIG. 10: Test of the von Schweidler law, Eq. (|A1|) . It is well verified roughly over two decades, for 
T < 250 K. Data for higher T depart from this law at short times. Symbols are as in Fig.^ Dashed 
line is the fit to Eq. (jXTJ of the data for T = 230 K over the fitting range 0.01 < A - C\ < 0.50 
with the result b = 0.6 ± 0.1. 
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